{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to ITK Segmentation in SimpleITK Notebooks\n",
    "\n",
    "<b>Goal</b>: To become familiar with basic segmentation algorithms available in ITK, and interactively explore their parameter space.\n",
    "\n",
    "Image segmentation filters process an image to partition it into (hopefully) meaningful regions. The output is commonly an image of integers where each integer can represent an object. The value 0 is commonly used for the background, and 1 (sometimes 255) for a foreground object.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(SimpleITK)\n",
    "\n",
    "source(\"downloaddata.R\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by loading our data and looking at it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "img_T1 <- ReadImage(fetch_data(\"nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT1.nrrd\"))\n",
    "\n",
    "# To visualize the labels image in color (RGB) we need to rescale the intensities to [0-255]\n",
    "img_T1_255 <- Cast(RescaleIntensity(img_T1), \"sitkUInt8\")\n",
    "\n",
    "Show(img_T1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholding\n",
    "\n",
    "Thresholding is the most basic form of segmentation. It simply labels the pixels of an image based on the intensity range without considering geometry or connectivity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "seg <- img_T1>200\n",
    "\n",
    "#Overlay segmentation onto the original image using the default alpha blending value of 0.5.\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "# A slightly more targeted thresholding approach with the values of interest inside a range\n",
    "seg <- BinaryThreshold(img_T1, lowerThreshold=100, upperThreshold=400, \n",
    "                       insideValue=1, outsideValue=0)\n",
    "\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selecting thresholds\n",
    "SimpleITK has a number of histogram based methods for automatic threshold selection for a bimodal distribution. These include <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1HuangThresholdImageFilter.html\">Huang</a>, <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1MaximumEntropyThresholdImageFilter.html\">MaximumEntropy</a>, <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1TriangleThresholdImageFilter.html\">Triangle</a>, and the popular <a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1OtsuThresholdImageFilter.html\">Otsu's</a> method. These methods create a histogram then use a heuristic to determine the threshold value which separates the foreground from background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "otsu_filter <- OtsuThresholdImageFilter()\n",
    "otsu_filter$SetInsideValue(0)\n",
    "otsu_filter$SetOutsideValue(1)\n",
    "seg <- otsu_filter$Execute(img_T1)\n",
    "\n",
    "cat(otsu_filter$GetThreshold())\n",
    "\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Region Growing Segmentation\n",
    "\n",
    "The first step of improvement upon the naive thresholding is a class of algorithms called region growing. These include:\n",
    "<ul>\n",
    "  <li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ConnectedThresholdImageFilter.html\">ConnectedThreshold</a></li>\n",
    "  <li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ConfidenceConnectedImageFilter.html\">ConfidenceConnected</a></li>\n",
    "  <li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1VectorConfidenceConnectedImageFilter.html\">VectorConfidenceConnected</a></li>\n",
    "  <li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1NeighborhoodConnectedImageFilter.html\">NeighborhoodConnected</a></li>\n",
    "</ul>\n",
    "\n",
    "We used an external program, 3D Slicer, to determine that index (132,142,96) was a good seed for the left lateral ventricle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "seed <- c(132,142,96)\n",
    "\n",
    "# visualize the seed with the original image: create a \"segmentation\" image with the seed \n",
    "# dilated so that it is clearly visible when overlaid onto the original image\n",
    "seg <- Image(img_T1$GetSize(), \"sitkUInt8\")\n",
    "seg$CopyInformation(img_T1)\n",
    "seg$SetPixel(seed,1)\n",
    "seg <- BinaryDilate(seg, 3)\n",
    "\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now segment using the connected threshold functionality, pixels that are connected to the seed(s) and lie within a given range of values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "seg <- ConnectedThreshold(img_T1, seedList=list(seed), lower=100, upper=190)\n",
    "\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Improving upon this is the ConfidenceConnected filter, which uses the initial seed(s) or current segmentation to estimate the thresholds based on the intensity mean and standard deviation: $\\mu\\pm c\\sigma$. \n",
    "\n",
    "The constant $c$ from the formula above is the \"multiplier\" function parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "seg <- ConfidenceConnected(img_T1, seedList=list(seed),\n",
    "                           numberOfIterations=1,\n",
    "                           multiplier=2.5,\n",
    "                           initialNeighborhoodRadius=1,\n",
    "                           replaceValue=1)\n",
    "\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have a multi-channel image we can possibly improve on the results above by using more channels. In our case we can use both T1 and T2 images to perform the segmentation. This is a multi-dimensional version of the previous filter with the \"multiplier\" used in conjunction with the Mahalanobis distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "img_T2 <- ReadImage(fetch_data(\"nac-hncma-atlas2013-Slicer4Version/Data/A1_grayT2.nrrd\"))\n",
    "\n",
    "img_multi <- Compose(img_T1, img_T2)\n",
    "seg <- VectorConfidenceConnected(img_multi, seedList=list(seed),\n",
    "                                 numberOfIterations=1,\n",
    "                                 multiplier=2.5,\n",
    "                                 initialNeighborhoodRadius=1)\n",
    "\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fast Marching Segmentation\n",
    "\n",
    "The FastMarchingImageFilter implements a fast marching solution to a simple level set evolution problem (eikonal equation). In this example, the speed term used in the differential equation is provided in the form of an image. The speed image is based on the gradient magnitude and mapped with the bounded reciprocal $1/(1+x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "seed <- c(132,142,96)\n",
    "feature_img <- GradientMagnitudeRecursiveGaussian(img_T1, sigma=.5)\n",
    "speed_img <- BoundedReciprocal(feature_img) # compute 1/(1+x) for each pixel\n",
    "\n",
    "Show(speed_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the FastMarchingImageFilter is a <b>time-crossing map</b> that indicates, for each pixel, how much time it would take for the front to arrive at the pixel location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "fm_filter <- FastMarchingBaseImageFilter()\n",
    "fm_filter$SetTrialPoints(list(seed))\n",
    "fm_filter$SetStoppingValue(1000)\n",
    "fm_img <- fm_filter$Execute(speed_img)\n",
    "\n",
    "Show(Threshold(fm_img,\n",
    "               lower=0.0,\n",
    "               upper=fm_filter$GetStoppingValue(),\n",
    "               outsideValue=fm_filter$GetStoppingValue()+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell you need to select the arrival time to obtain the desired segmentation (we start at 500 which is an over-segmentation, zero is the lower bound - rerun this cell to search for the desired threshold - binary search comes to mind)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "arrival_time <- 500\n",
    "seg <- fm_img<arrival_time\n",
    "\n",
    "Show(LabelOverlay(img_T1_255, seg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Level-Set Segmentation\n",
    "\n",
    "There are a variety of level-set based segmentation filters available in SimpleITK:\n",
    "<ul>\n",
    "<li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1GeodesicActiveContourLevelSetImageFilter.html\">GeodesicActiveContour</a></li>\n",
    "<li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ShapeDetectionLevelSetImageFilter.html\">ShapeDetection</a></li>\n",
    "<li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ThresholdSegmentationLevelSetImageFilter.html\">ThresholdSegmentation</a></li>\n",
    "<li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1LaplacianSegmentationLevelSetImageFilter.html\">LaplacianSegmentation</a></li>\n",
    "<li><a href=\"https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ScalarChanAndVeseDenseLevelSetImageFilter.html\">ScalarChanAndVese</a></li>\n",
    "</ul>\n",
    "\n",
    "First we create a label image from our seed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seed <- c(132,142,96)\n",
    "\n",
    "seg <- Image(img_T1$GetSize(), \"sitkUInt8\")\n",
    "seg$CopyInformation(img_T1)\n",
    "seg$SetPixel(seed, 1)\n",
    "seg <- BinaryDilate(seg, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the seed to estimate a reasonable threshold range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stats <- LabelStatisticsImageFilter()\n",
    "# ITK and therefore SimpleITK return an empty image for this execute - we'll just ignore it\n",
    "dev_null <- stats$Execute(img_T1, seg)\n",
    "\n",
    "factor <- 3.5\n",
    "lower_threshold <- stats$GetMean(label=1)-factor*stats$GetSigma(label=1)\n",
    "upper_threshold <- stats$GetMean(label=1)+factor*stats$GetSigma(label=1)\n",
    "cat(sprintf(\"lower threshold: %f\\nupper threshold %f\\n\",lower_threshold,upper_threshold))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "init_ls <- SignedMaurerDistanceMap(seg, TRUE, TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lsFilter <- ThresholdSegmentationLevelSetImageFilter()\n",
    "lsFilter$SetLowerThreshold(lower_threshold)\n",
    "lsFilter$SetUpperThreshold(upper_threshold)\n",
    "lsFilter$SetMaximumRMSError(0.02)\n",
    "lsFilter$SetNumberOfIterations(1000)\n",
    "lsFilter$SetCurvatureScaling(.5)\n",
    "lsFilter$SetPropagationScaling(1)\n",
    "lsFilter$ReverseExpansionDirectionOn()\n",
    "ls <- lsFilter$Execute(init_ls, Cast(img_T1, \"sitkFloat32\"))\n",
    "print(lsFilter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "Show(LabelOverlay(img_T1_255, ls>0))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
